ct = data3d('../data/lung/ct/00', 'short')
step = float_param(0.002, 0.0005, 0.01, 'S')

# rpl params
colormap = data1d_rgba('jet.png', 'uchar')
table = data1d('density-to-rpl-table.png', 'unsigned char')
rpl_threshold = float_param(4, 0, 4, 'T')
rpl_max = float_param(0.6, 0, 4, 'M')

# transfer function params
tf = data1d_rgba('bone.png', 'uchar')
tf_pos = float_param(0, -200, 200, 'P')
tf_width = float_param(100, 10, 400, 'W')
tf_opacity = float_param(1, 0, 1, 'O')

# lighting params
ambient = float_param(0.1, 0, 1, 'A')
lightpos = float3(-3, 3, 3)

cut = float_param(-2, -2, 2, 'C')

    
C = normalize(E - S)
result = float3(0)
see_through = 1.0
steps = length(E - S) / step
rpl = 0.0
t = rand() * step

for i in range(steps):
    P = (1-t) * S + t * E
    t += step
    
    # update radiological path length
    density = linear_query_3d(ct, P) * 32768
    rel_density = 0.21 + (1.78 - 0.21) * linear_query_1d(table,
        (density - (-766)) / (1173 - (-766)) * 2 - 1)
    rpl += rel_density * step

    # apply transfer function
    density = cubic_query_3d_cut(ct, P, C, cut) * 32768
    tf_query = (density - tf_pos) / tf_width
    if tf_query < 0: continue
    rgba = linear_query_1d_rgba(tf, tf_query*2 - 1)
    color = rgba.xyz
    
    # apply colormapped radiological path length
    if rpl < rpl_threshold:
        color = linear_query_1d_rgba(colormap, rpl / rpl_max * 2 - 1).xyz
    
    # compute lighting
    N = -normalize(cubic_gradient_3d_cut(ct, P, C, cut))
    L = normalize(lightpos - P)
    lit = phong(L, N, -C, color, 0.5, 50, ambient)

    # accumulate
    result += see_through * rgba.w * tf_opacity * lit
    see_through *= 1 - rgba.w * tf_opacity
    if see_through < 0.01: break

return result

